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A system of equations resulting from an approximation of the equation of motion of Green functions 
for correlated electron systems is usually solved using Matsubara technique. In this work we propose 
an alternative method which works entirely along the real frequency axis. Using the example of the 
attractive Hubbard model studied in the T-matrix approximation both self-consistently and non-self- 
consistently we demonstrate how powerful such a treatment is especially when dynamic quantities 
are calculated. 
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I. INTRODUCTION 

When investigating systems in thermodynamic equilib- 
rium with time independent Hamiltonians it was a great 
success to solve and discuss correlation functions in the 
transformed frequency space. These functions are ana- 
lytic functions in the complex plane except for a branch 
cut along the real axis. This knowledge was first used for 
the T=0 Green function technique e.g. pL§|- However, 
especially for finite temperatures, the theory of complex 
differentiable functions had led to the development of a 
very powerful method, the Matsubara technique || . 

For the derivation of the equations of motion an imag- 
inary axis formulation using Matsubara frequencies was 
used very early. The numerical solution of such equations 
however was first achieved along the real axis. An exam- 
ple for this is the solution of the Eliashberg equations by 
Schrieffer et al. §,§. 

Later it was discovered that a numerical solution was 
far more efficient when an imaginary axis technique is 
used. In this way static quantities like e.g. expectation 
values for occupation numbers (n) or double occupancy 
(n-|-nj.) and their temperature dependence and some ther- 
modynamic quantities can be calculated successfully. An 
example is the calculation of the the superconducting 
critical temperature |6|j^] and the temperature depen- 
dence of the critical magnetic field and specific heat from 
the Eliashberg equations || . 

A major shortcoming of such an imaginary axis treat- 
ment is that dynamical quantities like, for example, the 
electron density of states (along the real frequency axis) 
are difficult to obtain. Usually Pade approximants jjj 
or maximum entropy |i"c| techniques are used to obtain 
dynamical quantities. The fact that we need these com- 
plicated methods shows how difficult it is to calculate 
the values of a function at n R points along the real axis 
when n 1 points along the imaginary axis are known. The 
simplest way to illustrate this is the following: The fur- 
ther away the n 1 points are from the n R points the more 
the system of equations which has to be solved in or- 
der to obtain dynamical quantities tends to be singular 



and therefore tiny inaccuracies of the functions on the 
n 7 points, the Matsubara frequencies, can lead to much 
bigger inaccuracies along the real axis, at points n R . 

In this paper we go in some sense back to the old real 
axis technique and propose a treatment in which all quan- 
tities are directly calculated along the real axis. Since we 
do not calculate the functions at the Matsubara frequen- 
cies we do not need a mapping onto the real axis. As 
an example we use the T-matrix calculation using the 
ladder diagrams in the particle-particle channel for the 
attractive Hubbard model. This approximation is valid 
in the low density regime and is of particular interest 
since it might be able to model some aspects of the short 
coherence length pairs observed in the high-T c supercon- 
ductors 0HU|. 



II. SPECTRAL REPRESENTATIONS 

In order to motivate our method we give here a short 
overview of how different treatments of Green functions 
can be connected by using a spectral representation and 
by using the whole frequency plane. 

A one particle correlation function of operators C and 
B which are either both Bosonic or Fcrmionic operators 
can be written as a function of temporal and spatial co- 
ordinates: 

G k (x h t, Xj ,i/) = -i{T F/B tf{x u t)B{ Xj ,t')) (1) 

where Tpm is the time ordering operator for Fermions 
or Bosons and 



(...) = Z- x Tr(e- m '...) 



(2) 



is the thermal expectation value. Z is the partition func- 
tion in the grand canonical ensemble and H = H — fj,N is 
the Hamiltonian in this ensemble. If we restrict ourselves 
to systems in thermal equilibrium with periodic bound- 
ary conditions and time independent Hamiltonians we 
can apply a Fourier transform fTjfl: 
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G k (k,uj) = J°° d(t-t>)± J d( Xi - Xj -) 
G\t-t\ Xl 



Xj e 



iu:(t — t') ik(xi 



(3) 



The definition (|j]) and equation (||) are valid at all tem- 
peratures. 

The function G fe (k, uj) is related to a function which is 
analytic in the whole complex cj-plane with the exception 
of a branch cut along the real axis. Therefore there exists 
a spectral representation and G k (k, uj) can be rewritten 
as a function of a spectral function J(k, uj) which is a real 
function defined along the real uj axis. The Green func- 
tion G k (k,uj) is connected to its spectral representation 
in the following way: 



G k (k,uj) 



lim — 

<5->0+ 277 



ctUJ 



J(k,uJ)e 



— oc 
0tJ 



iS 



lj — UJ — i5 



(4) 



Here and in the following the upper sign corresponds to 
Fermions whereas the lower sign describes Bosons. For 
practical use we want to define a slightly different func- 
tion A(k, uj): 



A(k,uj) = y(k,uj)(e^±l) 



(5) 



Therefore equation (|]) can be rewritten: 
1 f°° 

G k (k,uj) = lim - / duJ 

-5-0+ 7T J_ ao 



A(k,ZJ)J^ A(k,U)- 



id 



id 



(6) 



For a non-interacting system and C,B e.g. Fermionic 
operators Ck the function ^ k A(k, uj) is identical to the 
density of states. For an interacting system J^k uj) 
will still be a density of states but A(k, uj) becomes de- 
pendent on the thermodynamic variables T,fi. 



A(k,tu) — >A T ^(k, 



(7) 



The usual T — Green functions for which the zero tem- 
perature diagram technique fl5|| is valid can be rewritten 
as a function of A T - ,1 (k, uj): 



G k (k,uj) 



lim lim ■ 

/3^oo g—,0+ 7T J_ r 



A T ^(k,uj) : 



dijj 



ss ±1 



UJ — UJ + i8 



± 



A T ^(k,uj) 



— uj — iS 



G R {k 1 uj) 
G A (k,uj) 



lim lim — 



lim lim — 

/3^oo 5-+0+ IT J_ r 



UJ — UJ + i5 

A T ^{k,uj) 

uj — uj — iS 



(8) 
(9) 
(10) 



where R, A denotes retarded and advanced Green func- 
tions. G R (k,u>) and G A (k, uj) are both branches of one 
function G(k, z) defined on the whole complex plane with 
the exception of the branch cut along the real axis where 
the poles of A T ' fJ, (k,ZJ) are located: 



G(k,z) = - 



duj 



A T ^(k,uj) 



(11) 



Also the thermal Green functions defined entirely along 
the imaginary frequency axis can be written as a function 
ofA T ^(k,w) 



7T / ™ lUJ n - UJ 



(12) 



where iuj r . 



(2 "+ 1} " , 2ip are the Matsubara frequen- 
cies for Fermions and Bosons respectively. The fact that 
the functions are only defined at certain periodic points 
means that the Fourier series of eq. ( fl2"| ) 

G 4 (k, ? r) = - e-^ T G\k,iuj n ) (13) 

71= — OO 

is periodic in imaginary time 

G t (k,iT + i/3) = TG t (k,ir) . (14) 

This means that the full knowledge of the function 
G(k, z) is either obtained by knowing (a) G(k, iuj n ) on the 
infinite but discrete points iuj n along the imaginary axis 
or by knowing (b) A T -^(k, uj) on a continuum along the 
real axis. The details of the usual Green function tech- 
nique are outlined in many textbooks e.g. fL6|]l5|| , Here 
we only want to highlight the connection of both meth- 
ods with the function A T ^(k : uj). The usual Matsubara 
technique determines the functions along the points iui n 
whereas we will discuss a method in this paper which 
calculates an approximation of A T ' fM (k, uj). 



III. NUMERICAL TECHNIQUE 

The solution of the equations in a certain approxima- 
tion for a correlated quantum system will generally be 
found numerically. In order to achieve self-consistency 
it is preferable to perform discrete sums over Matsubara 
frequencies than to calculate a function like A T,fJ, (k, uj). 
But to obtain dynamical quantities along the real axis, 
a difficult and somewhat uncontrolled analytic continua- 
tion will have to be performed. Therefore there have been 
attempts to solve such systems along the real axis (see for 
example [[l7| for the Eliashberg equations and [ p~8|JTg| ] for 
the self-consistent T-matrix equations). In these works 
some numerical integration was required along the real 
axis. 

However in this section we argue that it is possible to 
replace A T ' M (k, uj) by a series of (typically a few hundred) 
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S functions along the real axis. With this approximation 
all frequency integrations will turn into summations over 
a finite number of 6 functions and can therefore be done 
analytically. We use 



N„ 



1=1 



(15) 



where the amplitudes are af (which also depend on the 
thermodynamic variables) and the poles are located at 
the positions 6;. The Green function G(k, to) can now be 
expressed in this approximation as a sum of poles: 



G(k,KJ„) 



1 



iV mox 

E 

i=i 



(16) 



Using such a spectral representation our aim is to convert 
the usually complicated equations for G(k, iuj n ) into sets 
of equations for the amplitudes af only. 



A. Frequency grid 

An important point is how to choose the frequency 
points bi along the real axis. In a previous work |2(f | 
we let them fluctuate freely during the calculation but 
employed some approximations to restrict the number 
N m ax- In this work we keep them fixed relative to the 
chemical potential which leads to an efficient algorithm. 
It also turned out to be of importance to adopt the fre- 
quency grid to the problem (e.g. if the influence of a 
band edge is important the frequency points should be 
more dense around that band edge) and especially to the 
temperature. 

For a Fermionic system - where the energy range of 
±ksT around the chemical potential is of importance 



we choose for the example discussed in section IV 
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hi i2 = tanh 
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N„ 
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mm max 



N„ 



l 



(17) 



(18) 



and u) m in and w ma x are the minimum and maximum fre- 
quencies considered, respectively. The parameter a was 
adjusted in the way that the distance between two fre- 
quency points around the chemical potential is always 
smaller than fc#T. In this way the delta functions are ar- 
ranged with highest density at the chemical potential and 
far away from the chemical potential they get thinned 
out. Note that the index I does not need to have integer 
values - this will be used later on. 



B. Products of correlation functions 

Usually in such calculations we have to deal with prod- 
ucts of correlation functions which can be folded both 
in momentum and frequency space. The result of such 
products of one-particle correlation functions are general- 
ized susceptibilities. If we have two functions G 1 (k, iu> n ), 
G 2 (k, iu) n ), we often need to calculate the following prod- 
uct: 



xciwn)=~i;^£ 



(3^ N 

n v, 

G 1 (q, nu n ) G 2 (K - q, m m - iu3 n ) (19) 

where K is the total momentum of a pair of parti- 
cles. We denote in the following Bosonic Matsubara fre- 
quencies and q- vectors of Bosons with upper-case letters 
whereas for Fermionic systems we use lower-case letters. 
G l (k, iuj n ), for i = 1,2, has the approximate spectral 
representation: 



G*(q,iw„) 



JV ra „ 

E 



(20) 



We can directly evaluate the frequency summations and 
are left with sums over the coefficients. When inserting 
the spectral representation for G*(k, iuj n ) we get: 



2 K-q 

I 1 



j3 ^— ' iQ n - iu m — bi iu) m - bj 



EE 



fln-bj-bi \l + e^ 1 + , 



1 \ - >r-^ 



N a}*a 2K -<* 



ifl n - bj - bi 



tanh( Y +tanh ^)) (21) 



We now have determined a spectral representation of 
x(K, Z) which is valid in the whole complex plane. 
Nonetheless, for convenience we continue writing our 
function on the Matsubara frequencies. %(K, if2 n ) 
is however defined via a spectral representation on 
N max (N max + l)/2 frequency points along the real fre- 
quency axis which have to be folded back onto the N max 
points of our frequency grid. Therefore we sort the 
N max (N max + l)/2 points 



bf l =b j + b l , 



(22) 



check if they fit into each interval 6 p _!/ 2 < bj} < 

(the index of b in eq. ( |l7| ) can be non integer) and add 

their amplitudes 



cf t =a ; lq aj 2K - q (tanh' "" 



fa 
2 



+ tanh 



13k 



(23) 
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to the amplitude of the susceptibility. In this way we 
obtain a similar spectral representation for the suscepti- 
bility: 



x(K,ffl„)= J2 



P =i 



tilji bp 



(24) 



For the further calculation we just have to store the am- 
plitudes Cp . Note that if the G l (q,iuj n ),i = 1,2 had 
been Fermionic correlation functions, the amplitudes 
are - as opposed to the amplitudes a* k - not only posi- 
tive but change sign at the chemical potential due to the 
tanh function. The function x(K, iO, m ) will therefore be 
a Bosonic correlation function and iVl m are Bosonic Mat- 
subara frequencies. 



C. Vertex functions 



where t is the transfer integral between two neighbouring 
lattice sites < i,j > (we restrict ourselves to hyper-cubic 
lattices) and U < is the attractive interaction for two 
electrons on the same site of the lattice. In the case of 
U = the Hamiltonian can be diagonalized and gives the 
usual dispersion of a tight binding Hamiltonian e(k). 
The non- interacting Green function G^ (k, iu> n ) is ap- 
proximated by the pole bj which is nearest to e(k). 



G*°>(k,<. 



■Jn - bj 



b j-i/2 < e(k) < b j+1 



/2 



(28) 



(29) 



The ladder approximation p5[ takes into account all non- 
crossing scattering events of a pair in the particle-particle 
channel and becomes exact for all coupling strengths in 
the low density limit n — > 0. 



Vertex functions r(ki, k 2 , iw n \, ioJ n 2, ■ ■•) result from 
Bethe-Salpeter type equations and they are in general 
functionals of two particle prop agators similar to the 
one described in subsection III B , of a correlation U corr 



(which in general is a function of k and ui as well) and of 
one particle correlation functions. 

r(k 1 ,k 2 ,iw nl ,iw n2 ) = 

r{x(ki + k 2 iu> ni ),...} . (25) 

If they can be reduced to analytic functions in the com- 
plex plane they can be calculated by shifting the poles 
of the spectral representation at b n into the upper half 
plane by an amount 



(26) 



where the shift into the upper half plane depends on n 
and is smallest close to the chemical potential. With the 
help of such a vertex function a proper self-energy can 
be calculated in a way similar to the calculation of the 
susceptibility in subsection IIIB. In subsection IV A this 



is shown on the example of the attractive Hubbard model 
in the ladder approximation. 



IV. EXAMPLE: ATTRACTIVE HUBBARD 
MODEL IN THE LADDER APPROXIMATION 

The attractive Hubbard model shows superconductiv- 
ity in its weak coupling 3D limit. The range of intermedi- 
ate coupling and 2-3 dimensions is of particular interest 
for the high-T c cuprates [O] . The Hubbard Hamiltonian 




FIG. 1. schematic diagram of the non-self-consistent lad- 
der equations. The full thick lines represent full, interacting 
Green functions, the thin lines non interacting Green func- 
tions and the wavy lines are the Hubbard interaction U. 

This leads to the following system of equations |2l| ] 
which has to be solved either non self-consistently (Eqs. 
(pl|)-(p^)) which is diagrammatically shown in Fig. [I] 



GW(k,iw n ) = , — 

uv n - e(k) 

X (0) (K,^ m ) = -i£^E 



G (0) (q,iw„) G (0) (K-q,ffi 



u 



(k + q, iw n + iuJ m )G^ (q, tw m ) 
G" sc (k,*u,„) = G (o )(k)ia;n) -i_ E (o) (k)ia;n) 



(30) 

(31) 
(32) 



(33) 
(34) 



H = —t 



C i<7 C 3<r + U ^2 n iT n iU 



(27) 



<i,j>,a 



4 




FIG. 2. schematic diagram of the fully self-consistent lad- 
der equations. 



or self-consistently (Ens. (]3q)-(p8[)) which is diagram- 
matically shown in Fig. |2| 



f(K,iO n ) 
S(k,io; n ) 

G(k, iu n ) 



G(q, iw„) G(K - q, ifi TO - iu n ) (35) 
f/ 2 X (K,z^„) 



l-[/ X (K,iO r , 

m q 

T(k + q, iw„ + iw m )G(q, iu) m ) 
1 

G(°)(k,iw„)- 1 -I](k,^„) 



(36) 



(37) 
(38) 



In order to obtain self-consistency the Eqs. (|3q)-(|38|) 
have to be calculated iteratively until a stable self- 
consistent solution is obtained. 



A. Vertex function in the ladder approximation 

In the case of the ladder approximation the vertex 
function (Eqs. (^2[),(|36])) itself does not tend to zero as 

fl — > oo, instead limn >oo = U. However the function 

T(K, itt m ) can be decomposed into a function which ful- 
fills Kramers-Kronig relations T(K, i£l m ) plus a constant 
U ||. The function T(K,iQ m ) is U 2 times the two- 
particle propagator of a pair with total momentum K. 
In order to obtain a spectral representation for T(K, z) 
we have two possibilities, which we now discuss in turn. 



1. Complex evaluation 

In the following we show the evaluation of T(K, z) in 
the complex plane. Therefore we need an approximate 
expression for the complex function x(K, z) at the fre- 
quency points b m of our frequency grid of the real axis. 



x(KA 



n—1 



b n + iS n 



(39) 



with 



fin — 2(^ + 1/2 _ ^n-l/2) 

This has to be put into the equation for T 

U 2 X (K,b m ) 



r(K,6 m ) 



l-U X (K,b m ) 



(40) 



(41) 



The amplitudes g^ of T(K, b m ) at the points b m are then 
given by: 



25,, 



<-)„ 



■3(r(K,6 m )) 



(42) 



and we have obtained a spectral representation for 
r(K, b m ) which can be used for further calculation: 



N„ 



T(K,in m ) a J2 



<J 



K 



iQ m b n 



(43) 



2. Evaluation with partial fractions 

A second possibility to calculate the amplitudes is 
by rewriting Eq. ([y]) 



r(K,n) = 



(44) 
(45) 



and looking for the poles of the denominator. Since it is 
a polynomial of order N max this seems to be a hopeless 
procedure. However the roots are bracketed between the 
old poles b m and we know that the number of roots be- 
tween 6^ and b^ +1 is one except at zero frequency were 
it can not exceed two. Therefore this problem can be 
solved numerically. Having found the roots b^ we can 
evaluate the amplitudes by putting the solution into the 
numerator of Eq. (|i^). This procedure is just calculating 
partial fractions for eq. ([44[). An analogous procedure as 
for the susceptibility (see Eq. (|22|)) allows us to calcu- 
late the amplitudes g^ at the frequency points b m by 
adding all amplitudes for frequencies b^ in the interval 



-1/2 



-1/2- 



Both methods [V A 1 and IV A 2 work well. In 



we 



used method IV A 2 and for the calculations presented in 



this paper we use method IV A 1 . 



B. Self energy 

Having calculated the vertex function T(K, ifl m ) we 
can proceed and calculate the self-energy S'(k, iui n ). The 
' just reminds us that we do not calculate the full self- 
energy since we use T(K, iil m ) instead of T(K, iVL m ). The 
inclusion of the frequency independent part (U) of the 
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vertex function will lead to the Hartree term. For the 
frequency dependent self-energy, we have 

E'(k,w n ) = 

A? E ~R E r ( k + q ' iuJn + iuJm ^ G ( q ' iuJm ) 



N ^— ' /3 
q 



'eeIe 

q j,2 m 

, »„„ k+q q 

E £ - 91 "■' 



.9/ 



k+q 



iV 



3',! 



6/ 



1 



+ 



,/3bj 



1 



(46) 



The Bosonic distribution function is due to the Bosonic 
nature of a pair of Fermions described by r(K, zf2 m ). 
Again an analogous procedure as for the susceptibility 
(see Eq. (|22|)) yields the coefficients s k at the frequencies 
bj of our frequency grid. In this way we obtain a spectral 
representation for the self-energy 



N„ 



£'(k,iw„) = J2 



(47) 



a one-particle picture (the equation < N < 2 with 
N the total particle number is fulfilled for all possible 
values of fi and T). Already this is an improvement to 
the solution of the T-matrix equations given by Schmitt- 
Rink et al. |23) , as advocated by Serene Urn . 

It is also possible to go iteratively through the Eqs. 
(|35|)-(|38|). In this case the amplitudes will need an addi- 
tional index p for the number of iterations. 



,,K 



r YL,p 



(50) 



Now the Eqs. (g5|)-(gf 
next level of iteration. 



just map the amplitudes to the 



- a 



k,p+l 

'J 



K,p+1 



... (51) 

This procedure has to be repeated until a stable self- 
consistent solution is achieved. As a condition for self- 
consistency we used: 



N,; 



N„ 
k j=l 



1/2 



((« -' P ) 2 - (c 



< S (52) 



C. Calculation of the full Green function 

Knowing the self-energy S(k, iui n ) we can calculate the 
Green function for the one-particle propagator. 

G(k, iuj n ) = (iu n - e(k) +fi- S(k, iw„)) _1 
iu n — e(k) + 



JV m „ 

E 

3 = 1 



(48) 



Again we have two alternatives to calculate a spectral 
represe ntation for G(k,iuj n ). Analogous to the case 
IV A 1 when calculating the vertex function from the sus- 
ceptibility we get for the amplitudes of the one-particle 
Green function, 



2<L 



were S was typically chosen around 10 



V. NUMERICAL DISCUSSION OF SPECTRAL 
QUANTITIES 

In the following we will present several dynamical 
quantities which were obtained by using a numerical so- 
lution of the equations of the T-matrix in the ladder ap- 
proximation for the attractive Hubbard model on a 2D 
square lattice. 



•3(G(k,6 m )) 
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N„ 



^3 U m -e(k) +M - J2 



3 = 1 



i5i 



(49) 



where S n is defined in Eq. (|40|). Due to the approxima- 
tions made the sum rules might not always be fulfilled. 
Therefore it is necessary to perform a sum rule check and 
correct (if needed) the amplitudes a}^. 



D. Self-consistent vs. non-self-consistent 



Going once through the procedure Eqs. (|3l|)-(|34|) we 
obtain a non-self-consistent result which is conserving in 
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FIG. 3. For a 8x8 cubic 2D lattice, k B T = 0.55[t], 
U — — 4[t], n — — 2[t] and n w 0.7 the density of states is plot- 
ted, uj — is the position of the chemical potential. In (a) the 
full line corresponds to results obtained with Matsubara tech- 
nique and Pade approximates and the dashed line is the result 
of our spectral representation technique for Nmax = 300. In 
(b) we compare different numbers of N ma x = 100, 300, 500. 

In Figs. | and | we compare the density of states for 
a non-self-consistent calculation ( Eqs. (j3l|)-(|34|)) on an 
8x8 lattice. The temperature was chosen to 0.55[t] and 
the attractive interaction U = — 4[t] which is half the 
bandwidth of the non-interacting system. The chemical 
potential was fixed to fi = — 2[t] which led to an electron 
density of n ~ 0.7 at this temperature. 
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FIG. 4. For the same parameters as in Fig. ^ the k depen- 
dent spectral function along the (1, 1) direction is plotted, (a) 
shows the result from the Matsubara technique whereas (b-d) 
contain the result for N max = 100, 300, 500 respectively. 

For the spectral representation technique we used a 
frequency range of (see Eq. (|l8|)) u> m i n = — 24 [t] < lo < 
^max — 24 [t], the parameter a in Eq. (O) was cho- 
sen to be a = 2 and we discuss the effect of a different 
number of (^-functions N max . In Fig. ||a we compare 
the density of states as obtained with Matsubara tech- 
nique and a numerically exact analytic continuation onto 
the real axis |l^] which is possible only for the non-self- 
consistent calculation with a calculation for N max — 300. 
At to — fi = 2[t] there is a remnant of the logarithmic 
singularity which occurs in the middle of the band of a 
non-interacting 2D system. Below that, around uj = 
and below the chemical potential clear correlation effects 
can be seen which lead to additional states at u> < 0. In 
Fig. [SpD we compare for the same parameters different 
numbers N max of <5-peaks. For N max = 100 the correla- 
tion effects around the chemical potential are not clearly 
visible whereas for N max = 300 they are clearly present. 
Increasing N max up to 500 does not alter the picture. 

In Fig. |] we calculate k dependent quantities along 
the (1, 1) direction. The results from the Matsubara 
technique show that there is a strong incoherent broad- 
ening of the former quasiparticle peak around kp and 
for k < kp due to correlations. Along the diagonal the 
Fermi wave-vector is bracketed by (7r/4,7r/4) < kp < 
(7r/2,7r/2) Already the calculation with N max = 100 in 
Fig. ||b resolves the incoherent broadening but does fail 
to give further details which are clearly visible in the cal- 
culation for Fig. [|c for N max = 300 and in Fig. [|d for 

N max = 500. 

In order to demonstrate the strength of our method 
we discuss in the following (Figs. || - |l^) some aspects of 
the temperature dependence of the correlation functions 
obtained for a 16x16 lattice with N rnax = 300, u m i n = 
-32 [t] and uj rnax = 32[t] (see Eq. @). The strength of 
the attractive interaction for these figures is U = — 8[i], 
which is equal to the bandwidth of the non-interacting 
system. The particle number was chosen to be n = 0.2 
(n = 1 would correspond to half filling) and the chemical 
potential was adjusted as a function of temperature in 
order to keep the particle number n constant. In Figs, pi- 
ll non-self-consistent results according to Eqs. (^l|)-(§J) 
are presented whereas in Figs. |^-[l2], results from a fully 
self-consistent calculation are presented. 




FIG. 5. For a 16x16 cubic 2D lattice, U = -8[t] and 
n — 0.2 the density of states is plotted for three different 
temperatures, the chemical potential fj, has been adjusted as 
a function of temperature to keep n constant, u = is the 
position of the chemical potential, (a) shows the density of 
states for the three temperatures fcsT = 4.0[t] (dot-dashed 
line), k B T = 2.0[t] (full line) and k B T = 0.8[t] (dashed line), 
(b) shows for the lowest temperature the k-dependent spectral 
function along the (1, 1) direction. The results were obtained 
with a non-self-consistent calculation. 



Fig. ||a shows the k-integrated density of states. With 
decreasing temperature a gap occurs. The density at 
higher temperatures results from the one-particle con- 
tinuum and the density for low energies results from 
pairs in a two-particle bound state. Fig. [s|b shows that 
for kgT = 0.8[t] also the k-dependent spectral func- 
tion consists of two parts especially around k = kp, 
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(7r/8,7r/8). 




FIG. 6. For the same parameters as in Fig. |H| the imag- 
inary part of the susceptibility for pairs of electrons with 
total momentum K — is shown as obtained from a 
non-self-consistent calculation. 

Fig. ^ shows the imaginary part of the susceptibility 
for a total momentum of the pair of K = which is - as 
it should be - the non-interacting density of states where 
the energy has been stretched by a factor of 2. It is simply 
shifted according to the shift of the chemical potential 
with temperature. Due to the low density in the example 
we have chosen the chemical potential does not enter in 
the one-particle continuum. This is the case for higher 
densities where the imaginary part of the susceptibility 
changes sign at zero. 
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FIG. 7. For the same parameters as in Fig. ^ the imaginary 
part of the vertex function T(K = 0, ft) for pairs of electrons 
with total momentum K = is shown as it results from a 
non-self-consistent calculation. 

Fig. |^ shows the imaginary part of the vertex func- 
tion T(K = 0, ft). The strong peak corresponds to a 
true bound state and with decreasing temperature, the 
chemical potential drifts towards the bound state indicat- 
ing Bose condensation of non-interacting pairs into the 
bound state. 
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FIG. 8. For the same parameters as in Fig. ^ the imaginary 
part of the k-averaged self-energy is shown as it results from 
a non-self-consistent calculation. 
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In Fig. H we show the k-averaged imaginary part of the splitting into two parts, 
self-energy S(fc, ui). It mainly shifts with the chemical po- 
tential but does not otherwise show a large temperature 

dependence. 



0.2 



0.1 



0.0 



-10 




0.1 



k B T = 0.8 [t] 

■ k B T = 2.0 [t] 
k B T = 4.0 [t] 



a 

o 

II 



o.o 




-5 



10 15 20 25 30 

n [t] 



0.4 



0.2 



3 

-Si 



0.0 



k = (0,0) 

k = (7T/8.7T/8) 

k = (7r/47r/4) 
k = (w/2.w/2) 



FIG. 10. For the same parameters as in Fig. g the imagi- 
• 5 nary part of the susceptibility for pairs of electrons with to- 
tal momentum K — is shown as it results from a fully 
self-consistent calculation. 
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FIG. 9. For the same parameters as in Fig. g the den- 
sity of states resulting from a fully self-consistent calculation 
is plotted in (a), while (b) shows the k-dependent spectral 
function. 
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In the following figures we show results from a self- 
consistent calculation which have to be compared with 
the results from the non-self-consistent calculation. Fig. 

shows the density of states. Even though it looks 
similar for kgT — 4[<] there is strong difference at lower 
temperatures. The gap is no longer present due to the 
self-consistent procedure. Also for the k-dependent spec- 
tral functions (for k B T = Q.8[t] in Fig. |b) there is no 



FIG. 11. For the same parameters as in Fig. g the imagi- 
nary part of the vertex function T(K = 0, ft) for pairs of elec- 
trons with total momentum K — is shown as it results from 
a fully self-consistent calculation. 

Fig. [l^ shows the imaginary part of the susceptibility 
It is no longer an image of the non-interacting density 
of states. Note that it becomes negative for ft < at 
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the lowest temperature considered here (ksT = 0.8[t]) 
which is barely visible from the plot and indicates that 
the chemical potential is now in the one particle con- 
tinuum. Also the imaginary part of the vertex function 
T(K, f2), which is plotted for K = in Fig. [ll], changes 
sign at = 0, although this seems to happen in the plot 
(Fig. [ll]) for fl < 0, which is just an artifact of the broad- 
ening of the (5-functions which had to be applied in order 
to plot a spectral quantity. The peak in QF(K = 0, fi) no 
longer corresponds to a bound state; it is a two particle 
resonance in the one-particle continuum. 



the current paper. Here we have described in detail the 
numerical method and discussed its applicability to solve 
different problems of correlated quantum systems. 
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FIG. 12. For the same parameters as in Fig. ^ the imagi- 
nary part of the k-averaged self-energy is shown as it results 
from a fully self-consistent calculation. 

Also the self-energy whose imaginary part is (averaged 
over k) plotted in Fig. |l2| is strongly altered due to self- 
consistency. Compared to the non-self-consistent part in 
Fig. H it is strongly decreased in magnitude and starts 
at low temperatures to develop a minimum at u> = 
indicating the appearance of Fermi liquid like properties. 

VI. CONCLUSION 

Using the example of the attractive Hubbard model we 
have evaluated the ladder diagrams of the T-matrix, and 
we have demonstrated that our numerical method, which 
works entirely along the real frequency axis, enables us to 
accurately calculate spectral properties. We should point 
out that we did not show results for the lowest tempera- 
tures we were able to reach. In fact we can decrease the 
temperature for the calculations of Figs. || - [l2| by an 
additional two orders of magnitude without reaching nu- 
merical instabilities. However these results and especially 
their physical interpretation are not the main subject of 
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